Arbuscular mycorrhizal fungal communities in French vineyard: a signature of terroir whatever the farming practices

Code
knitr::opts_chunk$set(
  message = FALSE
)
Code
#if (!require("devtools", quietly = TRUE)) {
#  install.packages("devtools")
#}
#devtools::install_github("adrientaudiere/MiscMetabar")
Code
library(MiscMetabar)
library(targets)

library(metacoder)
library(DT)
library(knitr)
library(ggplot2)
library(patchwork)
library(vegan)
library(FactoMineR)
library(factoextra)
library(DESeq2)
library(indicspecies)
library(sf)
library(adespatial)
library(ade4)
library(adegraphics)
library(kableExtra)
library(formattable)
library(microbiomeMarker)
library(ComplexUpset)
Code
d <- tar_read(d_vs_blast)

Metadata modification

Code
d@sam_data$rank <- gsub("Weeding", "Grassing", d@sam_data$rank)
d@sam_data$inter_rank <- gsub("Weeding", "Grassing", d@sam_data$inter_rank)
d@sam_data$soil_practice <- gsub("Weeding", "Grassing", d@sam_data$soil_practice)

d@sam_data$terroir <- gsub("Perrier", "Vergèze", d@sam_data$terroir)

Spores number

Code
d@sam_data$nb_spores <- rowMeans(cbind(as.numeric(d@sam_data$spores_rep1), as.numeric(d@sam_data$spores_rep2),as.numeric(d@sam_data$spores_rep3)),na.rm=TRUE)

Practices

Code
d@sam_data$organic <- ifelse(!d@sam_data$practice %in%
  c("Conventional", "Conversion"),
"Organic", "NonOrganic"
)

Paired samples (roots, spores from the same samples)

Code
d@sam_data$paired_name <-
  apply(as.matrix(d@sam_data[, c(5:16)]), 1, paste,
    collapse = "--"
  )
d@sam_data$paired_name[d@sam_data$terroir == "Vergèze"] <-
  gsub("R", "", d@sam_data$ref_mycea[d@sam_data$terroir == "Vergèze"])
Code
tib_sam_data <- as_tibble(d@sam_data) %>%
  mutate(across(starts_with("Myc_freq"), as.numeric)) %>%
  mutate(across(starts_with("Myc_intensity_colonization"), as.numeric)) %>%
  mutate(across(starts_with("Myc_intensity_root"), as.numeric)) %>%
  mutate(across(starts_with("Arbuscul_richness"), as.numeric)) %>%
  mutate(across(starts_with("Arbuscul_abundance"), as.numeric)) %>%
  mutate(across(starts_with("Vesicle_richness"), as.numeric)) %>%
  mutate(across(starts_with("Vesicle_abundance"), as.numeric)) %>%
  mutate(Myc_freq = rowMeans(select(., starts_with("Myc_freq")), na.rm = TRUE)) %>%
  mutate(Myc_intensity_colonization = rowMeans(select(
    ., starts_with("Myc_intensity_colonization")
  ), na.rm = TRUE)) %>%
  mutate(Myc_intensity_root = rowMeans(select(., starts_with(
    "Myc_intensity_root"
  )), na.rm = TRUE)) %>%
  mutate(Arbuscul_richness = rowMeans(select(., starts_with(
    "Arbuscul_richness"
  )), na.rm = TRUE)) %>%
  mutate(Arbuscul_abundance = rowMeans(select(., starts_with(
    "Arbuscul_abundance"
  )), na.rm = TRUE)) %>%
  mutate(Vesicle_richness = rowMeans(
    select(., starts_with(
      "Vesicle_richness"
    )),
    na.rm = TRUE
  )) %>%
  mutate(Vesicle_abundance = rowMeans(
    select(., starts_with(
      "Vesicle_abundance"
    )),
    na.rm = TRUE
  )) %>% tibble::column_to_rownames(var = "X")

sample_data(d) <- tib_sam_data

Exclude Mycelium

Code
d <- clean_pq(subset_samples(d, compartment != "Mycelium"))
d_vs <- clean_pq(subset_samples(tar_read(d_vs), compartment != "Mycelium"))

Summary of final dataset

Code
summary_plot_pq(d) +
  labs(title = "Summary of final dataset") +
  theme(plot.title = element_text(hjust = 0.5, size = 20, color = "#1e2b4c"))

Summary of bioinformatics pipeline

Global summary

Code
kable(tar_read(track_sequences_samples_clusters))
nb_sequences nb_clusters nb_samples
Paired sequences 8283716 28637 194
Paired sequences without chimera 7777304 9791 194
Paired sequences without chimera and longer than 200bp 7736170 7522 194
ASV denoising 7533009 7522 186
OTU after vsearch reclustering at 97% 7533009 1081 186
OTU after blast filter without reclustering 3819912 3793 182
OTU after blast filter with reclustering 3822304 219 182

SUBSETING AND FILTERING DATASET

Code
d_muco <- clean_pq(subset_taxa(d, Class_PR2 == "Mucoromycota"))
d_vs_blast <- clean_pq(subset_taxa(tar_read(d_vs_blast), Class_PR2 == "Mucoromycota"))
Code
for(i in c(19:46)){
  d_muco@sam_data[,i] <- as.numeric(unlist(d_muco@sam_data[,i]))
}
Code
d_roots <- clean_pq(subset_samples(d_muco, compartment == "Roots"))
d_spores <- clean_pq(subset_samples(d_muco, compartment == "Spores"))
Code
d_rs_merged_samples <- clean_pq(merge_samples2(
  d_muco,
  "paired_name",
  default_fun = 
    function(x){MiscMetabar::diff_fct_diff_class(x, na.rm = T)}
))
Code
d_rs <- clean_pq(subset_samples_pq(
    d_rs_merged_samples,
    sample_sums(d_rs_merged_samples) > 2000
  ))

TAXONOMICAL ANALYSES

Code
# Figure S8
multitax_bar_pq(d_vs, "Supergroup_PR2",
  "Subdivision_PR2", "Class_PR2",
  nb_seq = TRUE
) +
  ggtitle("Number of sequences (log10) including non-Mucoromycota")

Number of sequences

Code
# Figure S6
dir.create("output/krona")
Warning in dir.create("output/krona"): 'output/krona' existe déjà
Code
krona(
  d_vs,
  ranks = c(11:19),
  "output/krona/OTU_post_clustered_allEuk",
  name = "OTU_post_clustered_allEuk"
)
krona(
  d,
  ranks = c(11:19),
  "output/krona/OTU_filtOnMaarjam_AM",
  name = "OTU_filtOnMaarjam_AM"
)
krona(
  d_muco,
  ranks = c(11:19),
  "output/krona/OTU_filtOnMaarjamAndPr2_AM",
  name = "OTU_filtOnMaarjamAndPr2_AM"
)
merge_krona(
  c(
    "output/krona/OTU_post_clustered_allEuk",
    "output/krona/OTU_filtOnMaarjam_AM",
    "output/krona/OTU_filtOnMaarjamAndPr2_AM"
  ),
  "output/krona/Euk_to_AM_filtering_nb_seq.html"
)

Number of OTU

Code
# Figure S7
krona(
  d_vs,
  ranks = c(11:19),
  "output/krona/OTU_post_clustered_allEuk_Nbotu",
      name = "OTU_post_clustered_allEuk",
  nb_seq = F
)
krona(
  d,
  ranks = c(11:19),
  "output/krona/OTU_filtOnMaarjam_AM_Nbotu",
    name = "OTU_filtOnMaarjam_AM",
  nb_seq = F
)
krona(
  d_muco,
  ranks = c(11:19),
  "output/krona/OTU_filtOnMaarjamAndPr2_AM_Nbotu",
  name = "OTU_filtOnMaarjamAndPr2_AM",
  nb_seq = F
)
merge_krona(
  c(
    "output/krona/OTU_post_clustered_allEuk_Nbotu",
    "output/krona/OTU_filtOnMaarjam_AM_Nbotu",
    "output/krona/OTU_filtOnMaarjamAndPr2_AM_Nbotu"
  ),
  "output/krona/Euk_to_AM_filtering_Nbotu.html"
)

NUMBER OF SPORES ANALYSES

Code
psm_otu <- psmelt(as_binary_otu_table(d_spores)) |>
  filter(Abundance > 0) |>
  group_by(Sample) |>
  summarise(
    "Abundance" = sum(Abundance),
    "region" = region[1],
    "nb_spores" = nb_spores[1]
  )
Code
psm_res <- psmelt_samples_pq(d_spores) %>% 
  mutate(nb_spores_log10 = log10(nb_spores))

Spores count <-> Alpha diversity

Code
# Fig S2
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_spores, rngseed = 221013)) %>% mutate(nb_spores_log10 = log10(nb_spores))

ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Abundance, type = "non-parametric") +
  ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Hill_0, type = "non-parametric") +
  ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Hill_1, type = "non-parametric") +
  ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Hill_2, type = "non-parametric")

Spores count <-> Practice

Code
# Figure S5
ggstatsplot::ggbetweenstats(psm_res, practice, nb_spores_log10, type = "non-parametric") 

Spores count <-> Terroir

Code
# Figure S4
ggstatsplot::ggbetweenstats(
  mutate(psm_res,
    terroir = reorder(psm_res$terroir, psm_res$nb_spores)
  ),
  terroir,
  nb_spores_log10,
  type = "non-parametric",
  centrality.plotting = F,
  package = "ggthemes",
  palette = "stata_economist"
) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

Spores count <-> Mycorrhization rate

Code
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Myc_freq, type = "non-parametric") +
  ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Myc_intensity_colonization, type = "non-parametric") +
  ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Arbuscul_abundance, type = "non-parametric")

MYCORRHIZATION MEASURE

Mycorrhization rate <-> Alpha diversity

F = Myc_freq

Code
psm_res <- psmelt_samples_pq(d_rs)
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_rs, rngseed = 221013))

# Figure S3a
(ggstatsplot::ggscatterstats(psm_res, Myc_freq, Hill_0, type = "non-parametrique") +
  ggstatsplot::ggscatterstats(psm_res, Myc_freq, Hill_1, type = "non-parametrique") +
  ggstatsplot::ggscatterstats(psm_res, Myc_freq, Hill_2, type = "non-parametrique")) /
  (ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_freq, Hill_0, type = "non-parametrique") +
    ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_freq, Hill_1, type = "non-parametrique") +
    ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_freq, Hill_2, type = "non-parametrique"))

M = Myc_intensity_colonization

Code
# Figure S3b
(ggstatsplot::ggscatterstats(psm_res, Myc_intensity_colonization, Hill_0, type = "non-parametrique") +
  ggstatsplot::ggscatterstats(psm_res, Myc_intensity_colonization, Hill_1, type = "non-parametrique") +
  ggstatsplot::ggscatterstats(psm_res, Myc_intensity_colonization, Hill_2, type = "non-parametrique")) /
  (ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_intensity_colonization, Hill_0, type = "non-parametrique") +
    ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_intensity_colonization, Hill_1, type = "non-parametrique") +
    ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_intensity_colonization, Hill_2, type = "non-parametrique"))

A = Arbuscul_abundance

Code
# Figure S3c
(ggstatsplot::ggscatterstats(psm_res, Arbuscul_abundance, Hill_0, type = "non-parametrique") +
  ggstatsplot::ggscatterstats(psm_res, Arbuscul_abundance, Hill_1, type = "non-parametrique") +
  ggstatsplot::ggscatterstats(psm_res, Arbuscul_abundance, Hill_2, type = "non-parametrique")) /
  (ggstatsplot::ggscatterstats(psm_res_rarefied, Arbuscul_abundance, Hill_0, type = "non-parametrique") +
    ggstatsplot::ggscatterstats(psm_res_rarefied, Arbuscul_abundance, Hill_1, type = "non-parametrique") +
    ggstatsplot::ggscatterstats(psm_res_rarefied, Arbuscul_abundance, Hill_2, type = "non-parametrique"))

Mycorrhization rate <-> terroir

Code
ggstatsplot::ggbetweenstats(psm_res_rarefied, terroir, Myc_freq, type = "non-parametrique", centrality.plotting = F) +
  ggstatsplot::ggbetweenstats(psm_res_rarefied, terroir, Myc_intensity_colonization, type = "non-parametrique", centrality.plotting = F) +
  ggstatsplot::ggbetweenstats(psm_res_rarefied, terroir, Arbuscul_abundance, type = "non-parametrique", centrality.plotting = F)

Mycorrhization rate <-> practice

Code
ggstatsplot::ggbetweenstats(psm_res_rarefied, practice, Myc_freq, type = "non-parametrique", centrality.plotting = F) +
  ggstatsplot::ggbetweenstats(psm_res_rarefied, practice, Myc_intensity_colonization, type = "non-parametrique", centrality.plotting = F) +
  ggstatsplot::ggbetweenstats(psm_res_rarefied, practice, Arbuscul_abundance, type = "non-parametrique", centrality.plotting = F)

ECOLOGICAL ANALYSES : ALPHA DIVERSITY

Compartment

Code
psm_res <- psmelt_samples_pq(d_muco)
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_muco, rngseed = 221013))
Code
ggstatsplot::ggbetweenstats(psm_res_rarefied,
  compartment,
  Hill_0,
  type = "non-parametrique"
) +
  ggstatsplot::ggbetweenstats(psm_res_rarefied,
    compartment,
    Hill_1,
    type = "non-parametrique"
  ) +
  ggstatsplot::ggbetweenstats(psm_res_rarefied,
    compartment,
    Hill_2,
    type = "non-parametrique"
  )

Code
ggplot(psm_res_rarefied, aes(y=Hill_0, x=compartment, color=practice))+
  geom_point() + 
  geom_line(aes(group = paired_name)) +
  facet_wrap(~terroir)

Terroir

Code
psm_res <- psmelt_samples_pq(d_rs)
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_rs, rngseed = 221013))
Code
# Figure 4
 (ggstatsplot::ggbetweenstats(
    mutate(
      psm_res_rarefied,
      terroir = reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_0)
    ),
    terroir,
    Hill_0,
    type = "non-parametrique",
    centrality.plotting = F,
    package = "ggthemes",
    palette = "stata_economist",
    point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
    0.9, size = 3, stroke = 0, na.rm = TRUE),
    violin.args = list(width = 0, linewidth = 0)
  )+ 
  scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
                             "#e15968","#784b43","#c200ab","#00b3f0","#953726",
                             "#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_0))),levels(as.factor(psm_res_rarefied$terroir)))])  + coord_flip()) /
 (ggstatsplot::ggbetweenstats(
    mutate(
      psm_res_rarefied,
      terroir = reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_1)
    ),
    terroir,
    Hill_1,
    type = "non-parametrique",
    centrality.plotting = F,
    package = "ggthemes",
    palette = "stata_economist",
    point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
    0.9, size = 3, stroke = 0, na.rm = TRUE),
    violin.args = list(width = 0, linewidth = 0)
  )+ 
  scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
                             "#e15968","#784b43","#c200ab","#00b3f0","#953726",
                             "#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_1))),levels(as.factor(psm_res_rarefied$terroir)))])  + coord_flip()) /
 (ggstatsplot::ggbetweenstats(
    mutate(
      psm_res_rarefied,
      terroir = reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_2)
    ),
    terroir,
    Hill_2,
    type = "non-parametrique",
    centrality.plotting = F,
    package = "ggthemes",
    palette = "stata_economist",
    point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
    0.9, size = 3, stroke = 0, na.rm = TRUE),
    violin.args = list(width = 0, linewidth = 0)
  )+ 
  scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
                             "#e15968","#784b43","#c200ab","#00b3f0","#953726",
                             "#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_2))),levels(as.factor(psm_res_rarefied$terroir)))])  + coord_flip()) + plot_layout(axes = "collect")

Code
# Figure S10
 (ggstatsplot::ggbetweenstats(
    mutate(
      psm_res,
      terroir = reorder(psm_res$terroir, psm_res$Hill_0)
    ),
    terroir,
    Hill_0,
    type = "non-parametrique",
    centrality.plotting = F,
    package = "ggthemes",
    palette = "stata_economist",
    point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
    0.9, size = 3, stroke = 0, na.rm = TRUE),
    violin.args = list(width = 0, linewidth = 0)
  )+ 
  scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
                             "#e15968","#784b43","#c200ab","#00b3f0","#953726",
                             "#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res$terroir, psm_res$Hill_0))),levels(as.factor(psm_res$terroir)))])  + coord_flip()) /
 (ggstatsplot::ggbetweenstats(
    mutate(
      psm_res,
      terroir = reorder(psm_res$terroir, psm_res$Hill_1)
    ),
    terroir,
    Hill_1,
    type = "non-parametrique",
    centrality.plotting = F,
    package = "ggthemes",
    palette = "stata_economist",
    point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
    0.9, size = 3, stroke = 0, na.rm = TRUE),
    violin.args = list(width = 0, linewidth = 0)
  )+ 
  scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
                             "#e15968","#784b43","#c200ab","#00b3f0","#953726",
                             "#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res$terroir, psm_res$Hill_1))),levels(as.factor(psm_res$terroir)))])  + coord_flip()) /
 (ggstatsplot::ggbetweenstats(
    mutate(
      psm_res,
      terroir = reorder(psm_res$terroir, psm_res$Hill_2)
    ),
    terroir,
    Hill_2,
    type = "non-parametrique",
    centrality.plotting = F,
    package = "ggthemes",
    palette = "stata_economist",
    point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
    0.9, size = 3, stroke = 0, na.rm = TRUE),
    violin.args = list(width = 0, linewidth = 0)
  )+ 
  scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
                             "#e15968","#784b43","#c200ab","#00b3f0","#953726",
                             "#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res$terroir, psm_res$Hill_2))),levels(as.factor(psm_res$terroir)))])  + coord_flip()) + plot_layout(axes = "collect")

Code
# Fig_7a
p_accu_balanced_mod_terroir <- 
  MiscMetabar::accu_plot_balanced_modality(d_rs,
                                           "terroir", 
                                           nperm = 999, 
                                           step = 300,
                                           quantile_prob=0.9,
                                           sample.size = 10000,
                                           progress_bar = FALSE)
Warning in vegan::rarefy(as.matrix(unclass(x[i, ])), n, se = TRUE): most
observed count data have counts 1, but smallest count is 3
Warning in vegan::rarefy(as.matrix(unclass(x[i, ])), n, se = TRUE): most
observed count data have counts 1, but smallest count is 3
Warning in vegan::rarefy(as.matrix(unclass(x[i, ])), n, se = TRUE): most
observed count data have counts 1, but smallest count is 3
Warning in vegan::rarefy(as.matrix(unclass(x[i, ])), n, se = TRUE): most
observed count data have counts 1, but smallest count is 17
Warning in max(dff$xlab, na.rm = TRUE): aucun argument pour max ; -Inf est
renvoyé
Code
p_accu_balanced_mod_terroir +
  xlab("Number of sequences") + 
  ylab("Number of OTUs") + 
  ggrepel::geom_label_repel(data=summarise(group_by(p_accu_balanced_mod_terroir$data,factor) ,y=max(X1,na.rm = T)),aes(label=factor,y=y,x=max(p_accu_balanced_mod_terroir$data$x,na.rm = T)), max.overlaps = 1900, force = 20) + xlim(c(1,25000)) + theme(legend.position = "none" )
Warning: Removed 14 rows containing missing values or values outside the scale range
(`geom_line()`).

Practice

Code
# Figure 5
ggstatsplot::ggbetweenstats(psm_res_rarefied,
  practice,
  Hill_0,
  centrality.plotting = F,
  type = "non-parametrique"
) + ylab("Hill 0 (richness)") +
  ggstatsplot::ggbetweenstats(psm_res_rarefied,
    practice,
    Hill_1,
    centrality.plotting = F,
    type = "non-parametrique"
  ) + ylab("Hill 1 (~Shannon)") +
  ggstatsplot::ggbetweenstats(psm_res_rarefied,
    practice,
    Hill_2,
    centrality.plotting = F,
    type = "non-parametrique"
  ) + ylab("Hill 2 (~Simpson)") +  plot_layout(axes = "collect")

Code
# Figure S11
ggstatsplot::ggbetweenstats(psm_res,
  practice,
  Hill_0,
  centrality.plotting = F,
  type = "non-parametrique"
) +
  ggstatsplot::ggbetweenstats(psm_res,
    practice,
    Hill_1,
    centrality.plotting = F,
    type = "non-parametrique"
  ) +
  ggstatsplot::ggbetweenstats(psm_res,
    practice,
    Hill_2,
    centrality.plotting = F,
    type = "non-parametrique"
  )

Code
# Fig_7b
p_accu_balanced_mod <- 
  MiscMetabar::accu_plot_balanced_modality(d_rs,
                                           "practice", 
                                           nperm = 999, 
                                           step = 300,
                                           quantile_prob=0.90,
                                           progress_bar = FALSE)
Warning in max(dff$xlab, na.rm = TRUE): aucun argument pour max ; -Inf est
renvoyé
Code
p_accu_balanced_mod +
  xlab("Number of sequences") +
  ylab("Number of OTUs") + 
  ggrepel::geom_label_repel(data=summarise(group_by(p_accu_balanced_mod$data,factor) ,y=max(X1,na.rm = T)),aes(label=factor,y=y,x=max(p_accu_balanced_mod$data$x,na.rm = T)), max.overlaps = 1900, force = 20) + xlim(c(1,18000)) + theme(legend.position = "none" )
Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_line()`).

ZOOM ON LOW-DIVERSITY SAMPLES

Code
otu_hill <- vegan::renyi(d_rs@otu_table, scale = c(0, 1, 2), hill = TRUE)

low_div_samples <-
  rownames(otu_hill[otu_hill$`0` < 10 |
    otu_hill$`1` < 2 |
    otu_hill$`2` < 2, ])

low_div_samples
 [1] "2019--Languedoc--Cevennes--St Maurice de Cazevieille--Conventional--Conventional--4.247388--44.035228--Conventional--Herbicide--Scratching--Herbicide"                   
 [2] "PE13"                                                                                                                                                                    
 [3] "2023--Champagne--Côte des Blancs--Chetillons--Chetillons Desh sous rang_plein chimique--Conventional--4.034283--48.945738--Conventional--Herbicide--Herbicide--Herbicide"
 [4] "2023--Champagne--Côte des Blancs--Chetillons--Chetillons Ecorces--Conventional--4.034333--48.945724--Conventional--Herbicide--Grassing--Herbicide"                       
 [5] "2023--Champagne--Côte des Blancs--Cramant--BIONNES JEB--Conversion--4.018985--48.988723--Conversion--Scratching--Grassing--Grassing"                                     
 [6] "2023--Champagne--Côte des Blancs--Cramant--BATEAU AN39 JEB--Conversion--4.021354--48.995715--Conversion--Scratching--Grassing--Grassing"                                 
 [7] "2023--Champagne--Côte des Blancs--Cramant--TDB AM29 JEB--Conversion--4.016860--48.996831--Conversion--Scratching--Grassing--Grassing"                                    
 [8] "2023--Champagne--Côte des Blancs--Mireaux--JEB VOISIN 1 MIREAUX DESSOUS--Conventional--3.994335--48.984188--Conventional--Herbicide--Herbicide--Herbicide"               
 [9] "2023--Champagne--Côte des Blancs--Avize--AVIZE Desh chimique en plein--Conventional--4.014914--48.976412--Conventional--Herbicide--Herbicide--Herbicide"                 
[10] "2023--Champagne--Côte des Blancs--Avize--AVIZE Travail meca sous rang--Organic--4.014914--48.976351--Organic grassed--Scratching--Grassing--Grassing"                    
[11] "2020--Champagne--Côte des Bar--Les Riceys--Conventional--Conventional--4.392132--48.009899--Conventional--Herbicide--Scratching--Herbicide"                              
[12] "2020--Bordelais--Langoiran--Langoiran--Conventional chimique--Conventional---0.356503--44.713204--Conventional--Herbicide--Herbicide--Herbicide"                         
[13] "2020--Champagne--Vallee de la Marne--Vincelles--Conventional--Conventional--3.931795--49.077467--Conventional--Herbicide--Scratching--Herbicide"                         
[14] "2020--Camargue--Montcalm--Montcalm--Conventional--Conventional--4.320967--43.584075--Conventional--Scratching--Scratching--Scratching"                                   
Code
sum(otu_hill$`0` < 10)
[1] 11
Code
sum(otu_hill$`1` < 2)
[1] 8
Code
sum(otu_hill$`2` < 2)
[1] 11
Code
length(low_div_samples)
[1] 14
Code
low_div_samples_table <- tibble(cbind(as_tibble(d_rs@sam_data[low_div_samples, c("terroir", "practice", "rank", "inter_rank", "compartment")]), apply(otu_hill[low_div_samples, ], 2, round, digits = 2), sample_sums(d_rs)[low_div_samples]))
Warning in class(x) <- tibble_class: La class(x) est une chaîne multiple
("tbl_df", "tbl", ...) ; le résultat ne sera plus un objet S4
Code
low_div_samples_table$compartment <- ifelse(is.na(low_div_samples_table$compartment), "Spores + Roots", "Spores")

colnames(low_div_samples_table) <- c("Terroir", "Global practice", "Rank practice", "Inter rank practice", "Compartment", "Hill 0 (richness)", "Hill 1 (~Shannon)", "Hill 2 (~Simpson)", "nb_seq")

low_div_samples_table <- arrange(low_div_samples_table, as.numeric(`Hill 0 (richness)`))

low_div_samples_table <-
  rbind(
    low_div_samples_table,
    c(
      "",
      "",
      "",
      "MEAN",
      "(low -diversity samples)",
      round(mean(low_div_samples_table$`Hill 0 (richness)`), 2),
      round(mean(low_div_samples_table$`Hill 1 (~Shannon)`), 2),
      round(mean(low_div_samples_table$`Hill 2 (~Simpson)`), 2),
      round(mean(low_div_samples_table$nb_seq), 2)
    ),
    c(
      "",
      "",
      "",
      "MEAN",
      "(all samples)",
      round(mean(otu_hill$`0`), 2),
      round(mean(otu_hill$`1`), 2),
      round(mean(otu_hill$`2`), 2),
      round(mean(sample_sums(d_rs)), 2)
    ),
     c(
      "",
      "",
      "",
      "MAX",
      "(all samples)",
      round(max(otu_hill$`0`), 2),
      round(max(otu_hill$`1`), 2),
      round(max(otu_hill$`2`), 2),
      round(max(sample_sums(d_rs)), 2)
    )
  )

# Table 3
kbl(low_div_samples_table) |>
  kable_classic(full_width = F, html_font = "Cambria")
Terroir Global practice Rank practice Inter rank practice Compartment Hill 0 (richness) Hill 1 (~Shannon) Hill 2 (~Simpson) nb_seq
Côte des Blancs Conversion Scratching Grassing Spores + Roots 3 1.87 1.68 42775
Côte des Blancs Conventional Herbicide Grassing Spores + Roots 4 1.71 1.39 28481
Langoiran Conventional Herbicide Herbicide Spores + Roots 4 2.38 1.87 12643
Côte des Blancs Organic Scratching Grassing Spores + Roots 5 1.99 1.49 35728
Vergèze Organic Scratching Ploughing Spores 6 4.33 3.65 6768
Côte des Blancs Conventional Herbicide Herbicide Spores + Roots 6 1.34 1.12 41991
Côte des Blancs Conversion Scratching Grassing Spores + Roots 7 2.74 2.13 49313
Côte des Blancs Conventional Herbicide Herbicide Spores + Roots 8 1.68 1.26 46114
Côte des Blancs Conversion Scratching Grassing Spores + Roots 8 1.17 1.05 45845
Vallee de la Marne Conventional Herbicide Scratching Spores + Roots 9 2.7 1.78 11777
Montcalm Conventional Scratching Scratching Spores 9 3.45 2.94 24675
Côte des Blancs Conventional Herbicide Herbicide Spores + Roots 10 1.45 1.16 45017
Cevennes Conventional Herbicide Scratching Spores + Roots 12 1.91 1.35 93267
Côte des Bar Conventional Herbicide Scratching Spores + Roots 17 2.43 1.51 28971
MEAN (low -diversity samples) 7.71 2.22 1.74 36668.93
MEAN (all samples) 19.87 6.23 4.35 48209.53
MAX (all samples) 47 15.06 9.04 101104

ECOLOGICAL ANALYSES : BETA DIVERSITY

Compartment

Code
# Figure S9
p <- phyloseq::plot_ordination(d_muco,
  vegan::decorana(vegdist(as(otu_table(d_muco), "matrix"),
    method = "robust.aitchison"
  )),
  color = "region",
  shape = "compartment"
) +
  geom_point(size = 3) +
  stat_ellipse(inherit.aes = F, aes(x = DCA1, y = DCA2, linetype = compartment))
Warning in phyloseq::plot_ordination(d_muco,
vegan::decorana(vegdist(as(otu_table(d_muco), : `Ordination species/OTU/taxa
coordinate indices did not match `physeq` index names. Setting corresponding
coordinates to NULL.
Code
p + geom_line(data = p$data, aes(group = paired_name))

Practice and terroir

Code
p <-
  phyloseq::plot_ordination(d_rs,
    ordinate(d_rs,
      method = "NMDS",
      distance = "bray"
    ),
    color = "terroir",
    shape = "practice"
  ) + facet_wrap("region") + scale_shape_manual(values=c(25,23,24))
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2485558 
Run 1 stress 0.2596844 
Run 2 stress 0.2663072 
Run 3 stress 0.2567211 
Run 4 stress 0.256791 
Run 5 stress 0.2626177 
Run 6 stress 0.2636841 
Run 7 stress 0.2612171 
Run 8 stress 0.2652946 
Run 9 stress 0.262289 
Run 10 stress 0.2523773 
Run 11 stress 0.2583239 
Run 12 stress 0.2599351 
Run 13 stress 0.2562442 
Run 14 stress 0.409148 
Run 15 stress 0.2593146 
Run 16 stress 0.2658447 
Run 17 stress 0.2580753 
Run 18 stress 0.2626311 
Run 19 stress 0.261243 
Run 20 stress 0.2584539 
*** Best solution was not repeated -- monoMDS stopping criteria:
     2: no. of iterations >= maxit
    18: stress ratio > sratmax
Code
p_nmds <- p + geom_point(
  data = select(p$data, -c(region)),
  fill = "grey40",
  color = "grey20",
  size = 1,
  alpha = 0.5
) +
  geom_point(size = 3, aes(fill=terroir)) +
  ggtitle("NMDS on bray distance") + 
  scale_fill_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
                             "#e15968","#784b43","#c200ab","#00b3f0","#953726",
                             "#e09199","#379d30","#ed5019","#ff63b0")) +
  scale_color_manual(values=rep("grey20",14))
# Figure 9
p_nmds

Code
# Figure 6a
upset_pq(d_rs, "practice", taxa_fill = "Family",
         set_sizes=(
        upset_set_size()
         + geom_text(aes(label=after_stat(count)), hjust=-0.3, color="white", stat='count')
    ))

Code
ggvenn_pq(d_rs, "practice")

Code
# Figure 6b
upset_pq(d_rs, "terroir", taxa_fill = "Family", min_size = 2,  height_ratio = 0.6
 ,set_sizes=(
        upset_set_size()
         + geom_label(aes(label=after_stat(count)), hjust=-0.3, size=2.5,  stat='count')
    ))
Warning: Removed 61 rows containing non-finite outside the scale range
(`stat_count()`).

Global beta-div analysis

Spatial data

Code
library(geodist)
lat_lon <- as_tibble(d_rs@sam_data[, c("lat", "long")])
lat_lon$lat <- as.numeric(lat_lon$lat)
lat_lon$long <- as.numeric(lat_lon$long)
dist_spatial_meter <- as.dist(geodist(lat_lon, measure = "geodesic"),
  upper = FALSE
)
Code
lon_lat_rs <- d_rs@sam_data[, c("long", "lat")]
lon_lat_rs$long <- as.numeric(lon_lat_rs$long)
lon_lat_rs$lat <- as.numeric(lon_lat_rs$lat)
MEM <- dbmem(lon_lat_rs, MEM.autocor = "non-null")
Code
test_MEM <- moran.randtest(MEM, nrepet = 999)
test_MEM$pvalue_adjust <- p.adjust(test_MEM$pvalue, method = "BH")
Code
d_rs@sam_data$MEM_1 <- MEM[, 1]
d_rs@sam_data$MEM_2 <- MEM[, 2]
res_ado_spatial_robAit <- adonis_pq(
  d_rs,
  "MEM_1 + MEM_2 + practice + inter_rank + rank  + terroir",
  correction_for_sample_size = TRUE,
  dist_method = "robust.aitchison"
)

res_ado_spatial_robAit_rarefy <- adonis_pq(
  rarefy_even_depth(d_rs, rngseed = 626),
  "MEM_1 + MEM_2 + practice + inter_rank + rank  + terroir",
  correction_for_sample_size = FALSE,
  dist_method = "robust.aitchison"
)

res_ado_spatial_bray <- adonis_pq(
  d_rs,
  "MEM_1 + MEM_2 + practice + inter_rank + rank  + terroir",
  correction_for_sample_size = TRUE,
  dist_method = "bray"
)


res_ado_spatial_bray_rarefy <- adonis_pq(
  rarefy_even_depth(d_rs, rngseed = 626),
  "MEM_1 + MEM_2 + practice + inter_rank + rank  + terroir",
  correction_for_sample_size = FALSE,
  dist_method = "bray"
)

df <- data.frame(res_ado_spatial_bray)
df$names <-
  factor(
    rownames(df),
    levels = c(
      "sample_size",
      "MEM_1",
      "MEM_2",
      "practice",
      "inter_rank",
      "rank",
      "terroir",
      "Residual",
      "Total"
    )
  )
df <- df %>% filter(!names == "Total")
p1 <-
  ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
  geom_bar(stat = "identity") +
  scale_fill_viridis_c() +
  geom_text(
    label = case_when(
      df$Pr..F. < 0.001 ~ "***",
      df$Pr..F. < 0.01 ~ "**",
      df$Pr..F. < 0.05 ~ "*",
      df$Pr..F. > 0.05 ~ "ns"
    ),
    nudge_y = 0.01
  ) +
  ggtitle("Permanova on Bray")

df <- data.frame(res_ado_spatial_robAit)
df$names <-
  factor(
    rownames(df),
    levels = c(
      "MEM_1",
      "MEM_2",
      "practice",
      "inter_rank",
      "rank",
      "terroir",
      "Residual",
      "Total"
    )
  )
df <- df %>% filter(!names == "Total")
p2 <-
  ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
  geom_bar(stat = "identity") +
  scale_fill_viridis_c() +
  geom_text(
    label = case_when(
      df$Pr..F. < 0.001 ~ "***",
      df$Pr..F. < 0.01 ~ "**",
      df$Pr..F. < 0.05 ~ "*",
      df$Pr..F. > 0.05 ~ "ns"
    ),
    nudge_y = 0.01
  ) +
  ggtitle("Permanova on robust Aitchison")

df <- data.frame(res_ado_spatial_bray_rarefy)
df$names <-
  factor(
    rownames(df),
    levels = c(
      "MEM_1",
      "MEM_2",
      "practice",
      "inter_rank",
      "rank",
      "terroir",
      "Residual",
      "Total"
    )
  )
df <- df %>% filter(!names == "Total")

p3 <-
  ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
  geom_bar(stat = "identity") +
  scale_fill_viridis_c() +
  geom_text(
    label = case_when(
      df$Pr..F. < 0.001 ~ "***",
      df$Pr..F. < 0.01 ~ "**",
      df$Pr..F. < 0.05 ~ "*",
      df$Pr..F. > 0.05 ~ "ns"
    ),
    nudge_y = 0.01
  ) +
  ggtitle("Permanova on Bray after rarefaction")

df <- data.frame(res_ado_spatial_robAit_rarefy)
df$names <-
  factor(
    rownames(df),
    levels = c(
      "MEM_1",
      "MEM_2",
      "practice",
      "inter_rank",
      "rank",
      "terroir",
      "Residual",
      "Total"
    )
  )
df <- df %>% filter(!names == "Total")
p4 <-
  ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
  geom_bar(stat = "identity") +
  scale_fill_viridis_c() +
  geom_text(
    label = case_when(
      df$Pr..F. < 0.001 ~ "***",
      df$Pr..F. < 0.01 ~ "**",
      df$Pr..F. < 0.05 ~ "*",
      df$Pr..F. > 0.05 ~ "ns"
    ),
    nudge_y = 0.01
  ) +
  ggtitle("Permanova on robust Aitchison after rarefaction")

(p1 + ylim(0, 0.12) + p2 + ylim(0, 0.12)) /
  (p3 + ylim(0, 0.12) + p4 + ylim(0, 0.12))
Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).

Code
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), d_rs@sam_data$terroir))
Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value Pr(>F)
Groups    13 0.16217 0.012475  0.6941 0.7615
Residuals 61 1.09633 0.017973               
Code
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), d_rs@sam_data$practice))
Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value Pr(>F)
Groups     2 0.05732 0.028658  2.1729 0.1213
Residuals 72 0.94962 0.013189               
Code
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$terroir))
Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value Pr(>F)
Groups    13 0.16217 0.012475  0.6941 0.7615
Residuals 61 1.09633 0.017973               
Code
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$practice))
Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value Pr(>F)
Groups     2 0.05732 0.028658  2.1729 0.1213
Residuals 72 0.94962 0.013189               
Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$terroir))
Analysis of Variance Table

Response: Distances
          Df Sum Sq Mean Sq F value    Pr(>F)    
Groups    13  393.8 30.2920  9.5001 2.306e-10 ***
Residuals 61  194.5  3.1886                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$practice))
Analysis of Variance Table

Response: Distances
          Df Sum Sq Mean Sq F value Pr(>F)
Groups     2  39.37  19.685  1.4366 0.2445
Residuals 72 986.60  13.703               
Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$terroir))
Analysis of Variance Table

Response: Distances
          Df Sum Sq Mean Sq F value    Pr(>F)    
Groups    13  393.8 30.2920  9.5001 2.306e-10 ***
Residuals 61  194.5  3.1886                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$practice))
Analysis of Variance Table

Response: Distances
          Df Sum Sq Mean Sq F value Pr(>F)
Groups     2  39.37  19.685  1.4366 0.2445
Residuals 72 986.60  13.703               

Soil properties

Code
soil_prop <- as_tibble(d_rs@sam_data) |>
  select(paired_name, practice, organic, terroir, Coarse_sand:CEC)
soil_prop_num <- soil_prop |>
  select(-all_of(c("practice", "organic", "terroir", "CaO", "Total_sand", "Total_filt"))) |>
  tibble::column_to_rownames("paired_name") |>
  mutate(across(everything(), as.numeric)) |>
  tidyr::drop_na()
Code
pca_test_res <- PCAtest::PCAtest(soil_prop_num, varcorr = T, plot = F, counter = F)

Sampling bootstrap replicates... Please wait

Calculating confidence intervals of empirical statistics... Please wait

Sampling random permutations... Please wait

Comparing empirical statistics with their null distributions... Please wait

========================================================
Test of PCA significance: 17 variables, 44 observations
1000 bootstrap replicates, 1000 random permutations
========================================================

Empirical Psi = 55.8085, Max null Psi = 9.0868, Min null Psi = 3.6595, p-value = 0
Empirical Phi = 0.4530, Max null Phi = 0.1828, Min null Phi = 0.1160, p-value = 0

Empirical eigenvalue #1 = 7.28487, Max null eigenvalue = 2.88796, p-value = 0
Empirical eigenvalue #2 = 3.15651, Max null eigenvalue = 2.39696, p-value = 0
Empirical eigenvalue #3 = 2.59627, Max null eigenvalue = 2.09788, p-value = 0
Empirical eigenvalue #4 = 1.16046, Max null eigenvalue = 1.87236, p-value = 1
Empirical eigenvalue #5 = 0.93253, Max null eigenvalue = 1.59515, p-value = 1
Empirical eigenvalue #6 = 0.74202, Max null eigenvalue = 1.48039, p-value = 1
Empirical eigenvalue #7 = 0.39204, Max null eigenvalue = 1.32508, p-value = 1
Empirical eigenvalue #8 = 0.24218, Max null eigenvalue = 1.18435, p-value = 1
Empirical eigenvalue #9 = 0.19075, Max null eigenvalue = 1.0485, p-value = 1
Empirical eigenvalue #10 = 0.08962, Max null eigenvalue = 0.99123, p-value = 1
Empirical eigenvalue #11 = 0.07595, Max null eigenvalue = 0.85026, p-value = 1
Empirical eigenvalue #12 = 0.05693, Max null eigenvalue = 0.7427, p-value = 1
Empirical eigenvalue #13 = 0.0369, Max null eigenvalue = 0.65767, p-value = 1
Empirical eigenvalue #14 = 0.02507, Max null eigenvalue = 0.55583, p-value = 1
Empirical eigenvalue #15 = 0.01767, Max null eigenvalue = 0.48237, p-value = 1
Empirical eigenvalue #16 = 0.00018, Max null eigenvalue = 0.40241, p-value = 1
Empirical eigenvalue #17 = 5e-05, Max null eigenvalue = 0.35079, p-value = 1

PC 1 is significant and accounts for 42.9% (95%-CI:38.8-49.7) of the total variation
PC 2 is significant and accounts for 18.6% (95%-CI:16.6-25) of the total variation
PC 3 is significant and accounts for 15.3% (95%-CI:10.5-17.9) of the total variation

The first 3 PC axes are significant and account for 76.7% of the total variation

Variables 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 15, 16, and 17 have significant loadings on PC 1
Variables 4, 8, 10, 12, and 14 have significant loadings on PC 2
Variables 5, 6, and 7 have significant loadings on PC 3

Variables 8, 10, 13, 15, and 17 have significant correlations with PC 1
Variables 8, 10, 12, and 14 have significant correlations with PC 2
Variables , and  have significant correlations with PC 3
Code
pval <- c()
for (i in seq_len(length(pca_test_res$`Percentage of variation of empirical PC's`))) {
  obs <- pca_test_res$`Percentage of variation of empirical PC's`[i]
  null_model <- pca_test_res$`Percentage of variation of randomized data`[, i]

  pval[i] <- (sum(obs < null_model) + 1) / (1 + nrow(pca_test_res$`Percentage of variation of randomized data`))
}

pca_test_pval_adj <- p.adjust(pval, method = "BH")
pca_test_pval_adj
 [1] 0.005661006 0.005661006 0.005661006 1.000000000 1.000000000 1.000000000
 [7] 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
[13] 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
Code
res_pca <- dudi.pca(soil_prop_num, scannf = F, nf = 4)
Code
p_pca_variable <- fviz_pca_var(res_pca,
  col.var = "cos2",
  gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
  repel = TRUE
)

p_pca_variable_axe3_4 <- fviz_pca_var(res_pca,
  col.var = "cos2", axes = c(3, 4),
  gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
  repel = TRUE
)

res_pca_var <- get_pca_var(res_pca)

# Fig. 2
(p_pca_variable + p_pca_variable_axe3_4) / free((factoextra::fviz_eig(res_pca)) +
  (ggcorrplot::ggcorrplot(res_pca_var$cor)))

Code
practice <- unlist(lapply(strsplit(rownames(res_pca$tab), "--"), function(x) {
  x[5]
}))

terroir <- unlist(lapply(strsplit(rownames(res_pca$tab), "--"), function(x) {
  x[2]
}))

p_pca_ind_terroir <- fviz_pca_ind(res_pca,
  habillage = terroir,
  label = FALSE,
  repel = TRUE,
  mean.point = FALSE,
  addEllipses = TRUE,
  pointsize = 3,
  ellipse.type = "convex"
) +
  scale_shape_manual(values = rep(16, length(unique(terroir))))

p_pca_ind_terroir$data$terroir <- terroir
p_pca_ind_terroir$data$practice <- practice

p_pca_ind_terroir_practice <- p_pca_ind_terroir + geom_point(aes(x = x, y = y),
  shape = case_when(
    practice == "Organic" ~ 2,
    practice == "Conventional" ~ 6,
    practice == "Conversion" ~ 5
  ), size = 3.5,
  color = rgb(0, 0, 0, 0.7)
)


# value
kw1_terroir <- kruskal.test(res_pca$tab[, 1], terroir)
kw2_terroir <- kruskal.test(res_pca$tab[, 2], terroir)
kw3_terroir <- kruskal.test(res_pca$tab[, 3], terroir)

kw1_practice <- kruskal.test(res_pca$tab[, 1], practice)
kw2_practice <- kruskal.test(res_pca$tab[, 2], practice)
kw3_practice <- kruskal.test(res_pca$tab[, 3], practice)

# Figure 3
p_pca_ind_terroir_practice + 
  annotate("text", x = 2, y = 4.5, label = paste("Kruskal-Wallis (axe1 - terroir): pval=", round(as.numeric(kw1_terroir$p.value),4)), size = 3) + 
  annotate("text", x = 2, y = 4.25, label = paste("Kruskal-Wallis (axe2 - terroir): pval=", round(as.numeric(kw2_terroir$p.value),4)), size = 3) + 
  annotate("text", x = 2, y = 4, label = paste("Kruskal-Wallis (axe3 - terroir): pval=", round(as.numeric(kw3_terroir$p.value),4)), size = 3) + 
  annotate("text", x = -2, y = 4.5, label = paste("Kruskal-Wallis (axe1 - practice): pval=", round(as.numeric(kw1_practice$p.value),4)), size = 3) + 
  annotate("text", x = -2, y = 4.25, label = paste("Kruskal-Wallis (axe2 - practice): pval=", round(as.numeric(kw2_practice$p.value),4)), size = 3) + 
  annotate("text", x = -2, y = 4, label = paste("Kruskal-Wallis (axe3 - practice): pval=", round(as.numeric(kw3_practice$p.value),4)), size = 3)
Warning: Removed 21 rows containing missing values or values outside the scale range
(`geom_point()`).

Code
# Figure S1
ggplot(d_rs@sam_data, aes(x = practice, y = as.numeric(Cu), fill = practice)) +
  geom_violin() +
  geom_jitter()
Warning: Removed 23 rows containing non-finite outside the scale range
(`stat_ydensity()`).
Warning: Removed 23 rows containing missing values or values outside the scale range
(`geom_point()`).

Code
kruskal.test(as.numeric(d_rs@sam_data$Cu), d_rs@sam_data$practice)

    Kruskal-Wallis rank sum test

data:  as.numeric(d_rs@sam_data$Cu) and d_rs@sam_data$practice
Kruskal-Wallis chi-squared = 0.72149, df = 2, p-value = 0.6972
Code
d_with_pca <-
  clean_pq(subset_samples_pq(d_rs, d_rs@sam_data$paired_name %in% rownames(soil_prop_num)))

res_pca_ind <- get_pca_ind(res_pca)

if (!identical(match(rownames(res_pca_ind$coord), d_with_pca@sam_data$paired_name), sort(match(rownames(res_pca_ind$coord), d_with_pca@sam_data$paired_name)))) {
  stop("ERROR")
}

d_with_pca <- add_info_to_sam_data(d_with_pca, res_pca_ind$coord)
Code
dim_df <- as_tibble(d_with_pca@sam_data) |>
  select(c(starts_with("Dim"), "nb_seq", "nb_otu"))

cor_results <- correlation::correlation(dim_df)
cor_results %>%
  summary(redundant = TRUE) %>%
  plot()

Code
lon_lat_with_pca <- d_with_pca@sam_data[, c("long", "lat")]
lon_lat_with_pca$long <- as.numeric(lon_lat_with_pca$long)
lon_lat_with_pca$lat <- as.numeric(lon_lat_with_pca$lat)
MEM_with_pca <- dbmem(lon_lat_with_pca, MEM.autocor = "non-null")

d_with_pca@sam_data$MEM_1 <- MEM_with_pca[, 1]
d_with_pca@sam_data$MEM_2 <- MEM_with_pca[, 2]
Code
res_ado_spatial_soil_bray <-
  adonis_pq(
    d_with_pca,
    "MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
    correction_for_sample_size = TRUE
  )
res_ado_spatial_soil_robAit <-
  adonis_pq(
    d_with_pca,
    "MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
    correction_for_sample_size = TRUE,
    dist_method = "robust.aitchison"
  )

res_ado_spatial_soil_bray_rarefy <-
  adonis_pq(
    rarefy_even_depth(d_with_pca, rngseed = 626),
    "MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
    correction_for_sample_size = FALSE
  )
res_ado_spatial_soil_robAit_rarefy <-
  adonis_pq(
    rarefy_even_depth(d_with_pca, rngseed = 626),
    "MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
    correction_for_sample_size = FALSE,
    dist_method = "robust.aitchison"
  )

anova(betadisper(phyloseq::distance(d_rs, method = "bray"), d_rs@sam_data$practice))
Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value Pr(>F)
Groups     2 0.05732 0.028658  2.1729 0.1213
Residuals 72 0.94962 0.013189               
Code
anova(betadisper(phyloseq::distance(d_rs, method = "bray"), d_rs@sam_data$terroir))
Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value Pr(>F)
Groups    13 0.16217 0.012475  0.6941 0.7615
Residuals 61 1.09633 0.017973               
Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$practice))
Analysis of Variance Table

Response: Distances
          Df Sum Sq Mean Sq F value Pr(>F)
Groups     2  39.37  19.685  1.4366 0.2445
Residuals 72 986.60  13.703               
Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$terroir))
Analysis of Variance Table

Response: Distances
          Df Sum Sq Mean Sq F value    Pr(>F)    
Groups    13  393.8 30.2920  9.5001 2.306e-10 ***
Residuals 61  194.5  3.1886                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# TABLE 2
kbl(res_ado_spatial_bray) |>
  kable_classic(full_width = F, html_font = "Cambria")
Df SumOfSqs R2 F Pr(>F)
sample_size 1 1.0052204 0.0670420 7.319743 0.001
MEM_1 1 1.5103278 0.1007296 10.997797 0.001
MEM_2 1 0.9164972 0.0611247 6.673684 0.001
practice 2 0.4789872 0.0319455 1.743928 0.036
inter_rank 3 0.5442769 0.0362999 1.321092 0.128
rank 2 0.6046954 0.0403295 2.201614 0.004
terroir 13 2.9300449 0.1954160 1.641215 0.001
Residual 51 7.0038313 0.4671126 NA NA
Total 74 14.9938812 1.0000000 NA NA
Code
kbl(res_ado_spatial_soil_bray) |>
  kable_classic(full_width = F, html_font = "Cambria")
Df SumOfSqs R2 F Pr(>F)
sample_size 1 0.9095040 0.1022131 6.3231283 0.001
MEM_1 1 0.2091359 0.0235034 1.4539720 0.149
MEM_2 1 1.0896243 0.1224557 7.5753755 0.001
Dim.1 1 0.3845781 0.0432202 2.6736956 0.007
Dim.2 1 0.2413033 0.0271185 1.6776088 0.080
Dim.3 1 0.0934249 0.0104994 0.6495160 0.783
practice 2 0.2815596 0.0316426 0.9787409 0.508
inter_rank 2 0.1453520 0.0163352 0.5052642 0.957
rank 2 0.4576684 0.0514343 1.5909201 0.070
terroir 11 2.2092074 0.2482782 1.3962758 0.022
Residual 20 2.8767531 0.3232993 NA NA
Total 43 8.8981110 1.0000000 NA NA
Code
# TABLE S2
kbl(res_ado_spatial_bray_rarefy) |>
  kable_classic(full_width = F, html_font = "Cambria")
Df SumOfSqs R2 F Pr(>F)
MEM_1 1 1.5224735 0.1079197 11.273304 0.001
MEM_2 1 1.1172045 0.0791924 8.272450 0.001
practice 2 0.5880731 0.0416852 2.177222 0.011
inter_rank 3 0.3890128 0.0275750 0.960161 0.478
rank 2 0.5861573 0.0415495 2.170130 0.016
terroir 13 2.8818779 0.2042804 1.641473 0.003
Residual 52 7.0226637 0.4977978 NA NA
Total 74 14.1074628 1.0000000 NA NA
Code
# TABLE S3
kbl(res_ado_spatial_soil_bray_rarefy) |>
  kable_classic(full_width = F, html_font = "Cambria")
Df SumOfSqs R2 F Pr(>F)
MEM_1 1 0.3924837 0.0479272 2.8426837 0.014
MEM_2 1 1.0983511 0.1341225 7.9551458 0.001
Dim.1 1 0.4058745 0.0495624 2.9396712 0.008
Dim.2 1 0.2866462 0.0350031 2.0761235 0.054
Dim.3 1 0.1150841 0.0140532 0.8335323 0.581
practice 2 0.2932271 0.0358067 1.0618937 0.388
inter_rank 2 0.1439472 0.0175778 0.5212911 0.945
rank 2 0.4694158 0.0573216 1.6999441 0.061
terroir 11 2.0847081 0.2545690 1.3726490 0.056
Residual 21 2.8994281 0.3540566 NA NA
Total 43 8.1891660 1.0000000 NA NA
Code
# TABLE S4
kbl(res_ado_spatial_robAit_rarefy) |>
  kable_classic(full_width = F, html_font = "Cambria")
Df SumOfSqs R2 F Pr(>F)
MEM_1 1 118.30117 0.0313656 2.7789327 0.001
MEM_2 1 110.52944 0.0293050 2.5963722 0.002
practice 2 117.62397 0.0311860 1.3815125 0.043
inter_rank 3 93.00477 0.0246586 0.7282373 0.965
rank 2 154.47993 0.0409577 1.8143917 0.002
terroir 13 964.07334 0.2556078 1.7420303 0.001
Residual 52 2213.67750 0.5869192 NA NA
Total 74 3771.69011 1.0000000 NA NA
Code
# TABLE S5
kbl(res_ado_spatial_soil_robAit_rarefy) |>
  kable_classic(full_width = F, html_font = "Cambria")
Df SumOfSqs R2 F Pr(>F)
MEM_1 1 112.72637 0.0339610 1.6387773 0.008
MEM_2 1 112.45902 0.0338805 1.6348907 0.006
Dim.1 1 191.81334 0.0577875 2.7885165 0.001
Dim.2 1 93.30719 0.0281106 1.3564679 0.087
Dim.3 1 112.44158 0.0338752 1.6346370 0.014
practice 2 128.09782 0.0385919 0.9311211 0.678
inter_rank 2 98.34228 0.0296275 0.7148331 0.956
rank 2 158.00090 0.0476008 1.1484814 0.175
terroir 11 867.57636 0.2613741 1.1465934 0.057
Residual 21 1444.52439 0.4351909 NA NA
Total 43 3319.28925 1.0000000 NA NA

Variance partitionning

Code
# Figure 8a
res_varpart_rarefy <- var_par_rarperm_pq(
  physeq = d_with_pca,
  list_component = list(
    "Spatial" = c("MEM_1", "MEM_2"),
    "Soil" = c("Dim.1", "Dim.2", "Dim.3"),
    "Practice" = c("practice", "inter_rank", "rank"),
    "Terroir" = c("terroir")
  ),
  nperm = 99,
  dbrda_computation = TRUE,
  progress_bar = FALSE
)
plot_var_part_pq(res_varpart_rarefy, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)

Code
# Figure 8b
res_varpart_rarefy_wo_soil <- var_par_rarperm_pq(
  physeq = d_rs,
  list_component = list(
    "Spatial" = c("MEM_1", "MEM_2"),
    "Practice" = c("practice", "inter_rank", "rank"),
    "Terroir" = c("terroir")
  ),
  nperm = 99,
  dbrda_computation = TRUE,
  progress_bar = FALSE
)

plot_var_part_pq(res_varpart_rarefy_wo_soil, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)

Code
res_varpart_robAit_rarefy <- var_par_rarperm_pq(
  physeq = d_with_pca,
  list_component = list(
    "Spatial" = c("MEM_1", "MEM_2"),
    "Soil" = c("Dim.1", "Dim.2", "Dim.3"),
    "Practice" = c("practice", "inter_rank", "rank"),
    "Terroir" = c("terroir")
  ),
  dist_method = "robust.aitchison",
  nperm = 99,
  progress_bar = FALSE
)
plot_var_part_pq(res_varpart_robAit_rarefy, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)

Code
# Figure S12b
res_varpart_robAit_rarefy_wo_soil <- var_par_rarperm_pq(
  physeq = d_rs,
  list_component = list(
    "Spatial" = c("MEM_1", "MEM_2"),
    "Practice" = c("practice", "inter_rank", "rank"),
    "Terroir" = c("terroir")
  ),
    dist_method = "robust.aitchison",
  nperm = 99,
  dbrda_computation = TRUE,
  progress_bar = FALSE
)

plot_var_part_pq(res_varpart_robAit_rarefy_wo_soil, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)

Differential abundance analysis

Indicspecies

Terroir

Code
# Figure 10 
res_mpt_terroir <-
  multipatt(as.matrix(d_rs@otu_table),
    d_rs@sam_data$terroir,
    control = how(nperm = 9999),
    max.order = 3
  )

res_mpt_terroir_df <- res_mpt_terroir$sign
res_mpt_terroir_df$p.adj <- p.adjust(res_mpt_terroir_df$p.value, method = "BH")
res_mpt_terroir_df$OTU_names <- rownames(res_mpt_terroir_df)
res_mpt_terroir_df_signif <-
  res_mpt_terroir_df %>%
  filter(p.adj < 0.05) %>%
  tidyr::pivot_longer(cols = starts_with("s."))

tax_tab <- as.data.frame(d_rs@tax_table)
tax_tab$otu <- rownames(tax_tab)

res_mpt_terroir_df_signif_taxo <- left_join(res_mpt_terroir_df_signif,tax_tab, by = join_by("OTU_names" == "otu"))

ggplot(
  res_mpt_terroir_df_signif_taxo,
  aes(
    x = OTU_names,
    y = name,
    alpha = value,
    color = Genus
  )
) +
  geom_point(size = 6) +
  theme(axis.text.x = element_text(
    angle = 90,
    vjust = 0.5,
    hjust = 1
  ))+
  scale_alpha(range = c(0, 1))

Code
res_mpt_terroir_rg <-
  multipatt(as.matrix(rarefy_even_depth(d_rs, rngseed = 626)@otu_table),
    d_rs@sam_data$terroir,
    control = how(nperm = 9999),
    max.order = 3,
    func = "r.g"
  )

res_mpt_terroir_rg_df <- res_mpt_terroir_rg$sign
res_mpt_terroir_rg_df$p.adj <- p.adjust(res_mpt_terroir_rg_df$p.value, method = "BH")
res_mpt_terroir_rg_df$OTU_names <- rownames(res_mpt_terroir_rg_df)
res_mpt_terroir_rg_df_signif <-
  res_mpt_terroir_rg_df %>%
  filter(p.adj < 0.05) %>%
  tidyr::pivot_longer(cols = starts_with("s."))

ggplot(
  res_mpt_terroir_rg_df_signif,
  aes(
    x = OTU_names,
    y = name,
    alpha = value,
    color = stat
  )
) +
  geom_point(size = 4) +
  theme(axis.text.x = element_text(
    angle = 90,
    vjust = 0.5,
    hjust = 1
  ))

practice

Code
res_mpt_practice <-
  multipatt(as.matrix(d_rs@otu_table),
    d_rs@sam_data$practice,
    control = how(nperm = 9999),
    max.order = 3
  )

res_mpt_practice_df <- res_mpt_practice$sign
res_mpt_practice_df$p.adj <- p.adjust(res_mpt_practice_df$p.value, method = "BH")
res_mpt_practice_df$OTU_names <- rownames(res_mpt_practice_df)
res_mpt_practice_df_signif <-
  res_mpt_practice_df %>%
  filter(p.adj < 0.05) %>%
  tidyr::pivot_longer(cols = starts_with("s."))
res_mpt_practice_df_signif
# A tibble: 0 × 7
# ℹ 7 variables: index <int>, stat <dbl>, p.value <dbl>, p.adj <dbl>,
#   OTU_names <chr>, name <chr>, value <dbl>
Code
res_mpt_practice_rg <-
  multipatt(as.matrix(d_rs@otu_table),
    d_rs@sam_data$practice,
    control = how(nperm = 9999),
    max.order = 3,
    func = "r.g"
  )

res_mpt_practice_rg_df <- res_mpt_practice_rg$sign
res_mpt_practice_rg_df$p.adj <- p.adjust(res_mpt_practice_rg_df$p.value, method = "BH")
res_mpt_practice_rg_df$OTU_names <- rownames(res_mpt_practice_rg_df)
res_mpt_practice_rg_df_signif <-
  res_mpt_practice_rg_df %>%
  filter(p.adj < 0.05) %>%
  tidyr::pivot_longer(cols = starts_with("s."))

res_mpt_practice_rg_df_signif
# A tibble: 0 × 7
# ℹ 7 variables: index <int>, stat <dbl>, p.value <dbl>, p.adj <dbl>,
#   OTU_names <chr>, name <chr>, value <dbl>

Session information

Code
sessionInfo()
R version 4.3.3 (2024-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 11 (bullseye)

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.13.so;  LAPACK version 3.9.0

locale:
 [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
 [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
 [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Paris
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] geodist_0.0.8               ComplexUpset_1.3.3         
 [3] microbiomeMarker_1.8.0      formattable_0.2.1          
 [5] kableExtra_1.4.0            adegraphics_1.0-21         
 [7] ade4_1.7-22                 adespatial_0.3-23          
 [9] sf_1.0-16                   indicspecies_1.7.14        
[11] DESeq2_1.42.1               SummarizedExperiment_1.32.0
[13] Biobase_2.62.0              MatrixGenerics_1.14.0      
[15] matrixStats_1.3.0           GenomicRanges_1.54.1       
[17] GenomeInfoDb_1.38.8         IRanges_2.36.0             
[19] S4Vectors_0.40.2            BiocGenerics_0.48.1        
[21] factoextra_1.0.7            FactoMineR_2.11            
[23] vegan_2.6-6.1               lattice_0.22-6             
[25] permute_0.9-7               patchwork_1.2.0            
[27] knitr_1.46                  DT_0.33                    
[29] metacoder_0.3.7             targets_1.7.0              
[31] MiscMetabar_0.9.2           purrr_1.0.2                
[33] dplyr_1.1.4                 dada2_1.30.0               
[35] Rcpp_1.0.12                 ggplot2_3.5.1              
[37] phyloseq_1.46.0            

loaded via a namespace (and not attached):
  [1] igraph_2.0.3                    mia_1.10.0                     
  [3] Formula_1.2-5                   scater_1.30.1                  
  [5] rematch2_2.1.2                  zlibbioc_1.48.2                
  [7] tidyselect_1.2.1                bit_4.0.5                      
  [9] doParallel_1.0.17               BWStest_0.2.3                  
 [11] clue_0.3-65                     rjson_0.2.21                   
 [13] blob_1.2.4                      stringr_1.5.1                  
 [15] rngtools_1.5.2                  S4Arrays_1.2.1                 
 [17] parallel_4.3.3                  RNeXML_2.4.11                  
 [19] correlation_0.8.4               png_0.1-8                      
 [21] cli_3.6.3                       ggplotify_0.1.2                
 [23] CVXR_1.0-12                     bayestestR_0.13.2              
 [25] multtest_2.58.0                 MultiAssayExperiment_1.28.0    
 [27] bluster_1.12.0                  BiocNeighbors_1.20.2           
 [29] TreeSummarizedExperiment_2.10.0 adegenet_2.1.10                
 [31] mime_0.12                       evaluate_0.23                  
 [33] tidytree_0.4.6                  ComplexHeatmap_2.18.0          
 [35] ggh4x_0.2.8                     stringi_1.8.4                  
 [37] backports_1.4.1                 PMCMRplus_1.9.10               
 [39] lmerTest_3.1-3                  gsl_2.1-8                      
 [41] XML_3.99-0.16.1                 Exact_3.2                      
 [43] ggVennDiagram_1.5.2             httpuv_1.6.15                  
 [45] Wrench_1.20.0                   paletteer_1.6.0                
 [47] magrittr_2.0.3                  leaps_3.1                      
 [49] splines_4.3.3                   jpeg_0.1-10                    
 [51] doRNG_1.8.6                     wk_0.9.1                       
 [53] rootSolve_1.8.2.4               ggbeeswarm_0.7.2               
 [55] DBI_1.2.2                       withr_3.0.0                    
 [57] class_7.3-22                    systemfonts_1.0.6              
 [59] htmlwidgets_1.6.4               fs_1.6.4                       
 [61] SuppDists_1.1-9.7               SingleCellExperiment_1.24.0    
 [63] ggrepel_0.9.5                   labeling_0.4.3                 
 [65] SparseArray_1.2.4               cellranger_1.1.0               
 [67] flashClust_1.01-2               rncl_0.8.7                     
 [69] see_0.8.4                       lmom_3.0                       
 [71] effectsize_0.8.8                zoo_1.8-12                     
 [73] XVector_0.42.0                  decontam_1.22.0                
 [75] secretbase_0.5.0                foreach_1.5.2                  
 [77] fansi_1.0.6                     caTools_1.18.2                 
 [79] grid_4.3.3                      data.table_1.15.4              
 [81] ggtree_3.10.1                   rhdf5_2.46.1                   
 [83] irlba_2.3.5.1                   gridGraphics_0.5-1             
 [85] DescTools_0.99.54               base64url_1.4                  
 [87] lazyeval_0.2.2                  yaml_2.3.8                     
 [89] survival_3.6-4                  crayon_1.5.3                   
 [91] RColorBrewer_1.1-3              tidyr_1.3.1                    
 [93] later_1.3.2                     codetools_0.2-19               
 [95] base64enc_0.1-3                 GlobalOptions_0.1.2            
 [97] shape_1.4.6.1                   limma_3.58.1                   
 [99] estimability_1.5.1              Rsamtools_2.18.0               
[101] ggside_0.3.1                    foreign_0.8-86                 
[103] pkgconfig_2.0.3                 ggpubr_0.6.0                   
[105] xml2_1.3.6                      GenomicAlignments_1.38.2       
[107] aplot_0.2.2                     scatterplot3d_0.3-44           
[109] ape_5.8                         viridisLite_0.4.2              
[111] xtable_1.8-4                    interp_1.1-6                   
[113] car_3.1-2                       highr_0.10                     
[115] hwriter_1.3.2.1                 plyr_1.8.9                     
[117] httr_1.4.7                      rbibutils_2.2.16               
[119] tools_4.3.3                     broom_1.0.5                    
[121] beeswarm_0.4.0                  htmlTable_2.4.2                
[123] checkmate_2.3.1                 nlme_3.1-164                   
[125] PCAtest_0.0.1                   lme4_1.1-35.3                  
[127] digest_0.6.36                   numDeriv_2016.8-1.1            
[129] adephylo_1.1-16                 Matrix_1.6-5                   
[131] farver_2.1.2                    reshape2_1.4.4                 
[133] yulab.utils_0.1.4               viridis_0.6.5                  
[135] DirichletMultinomial_1.44.0     rpart_4.1.23                   
[137] glue_1.7.0                      cachem_1.0.8                   
[139] Hmisc_5.1-2                     generics_0.1.3                 
[141] Biostrings_2.70.3               classInt_0.4-10                
[143] mvtnorm_1.2-4                   spdep_1.3-3                    
[145] metagenomeSeq_1.43.0            statmod_1.5.0                  
[147] biomformat_1.30.0               ScaledMatrix_1.10.0            
[149] carData_3.0-5                   minqa_1.2.6                    
[151] ANCOMBC_2.4.0                   glmnet_4.1-8                   
[153] utf8_1.2.4                      seqinr_4.2-36                  
[155] gtools_3.9.5                    readxl_1.4.3                   
[157] datawizard_0.10.0               ggsignif_0.6.4                 
[159] gridExtra_2.3                   shiny_1.8.1.1                  
[161] GenomeInfoDbData_1.2.11         energy_1.7-11                  
[163] rhdf5filters_1.14.1             parameters_0.21.6              
[165] RCurl_1.98-1.14                 memoise_2.0.1                  
[167] rmarkdown_2.26                  scales_1.3.0                   
[169] gld_2.6.6                       svglite_2.1.3                  
[171] phylobase_0.8.12                rstudioapi_0.16.0              
[173] cluster_2.1.6                   hms_1.1.3                      
[175] ggcorrplot_0.1.4.1              munsell_0.5.1                  
[177] colorspace_2.1-0                rlang_1.1.4                    
[179] s2_1.1.6                        DelayedMatrixStats_1.24.0      
[181] sparseMatrixStats_1.14.0        circlize_0.4.16                
[183] scuttle_1.12.0                  mgcv_1.9-1                     
[185] xfun_0.43                       multcompView_0.1-10            
[187] coda_0.19-4.1                   e1071_1.7-14                   
[189] TH.data_1.1-2                   plotROC_2.3.1                  
[191] iterators_1.0.14                emmeans_1.10.1                 
[193] abind_1.4-5                     tibble_3.2.1                   
[195] treeio_1.26.0                   gmp_0.7-4                      
[197] Rhdf5lib_1.24.2                 DECIPHER_2.30.0                
[199] bitops_1.0-7                    Rdpack_2.6                     
[201] ps_1.7.6                        promises_1.3.0                 
[203] RSQLite_2.3.6                   sandwich_3.1-0                 
[205] DelayedArray_0.28.0             proxy_0.4-27                   
[207] Rmpfr_0.9-5                     compiler_4.3.3                 
[209] statsExpressions_1.5.4          prettyunits_1.2.0              
[211] boot_1.3-30                     beachmat_2.18.1                
[213] BiocSingular_1.18.0             units_0.8-5                    
[215] MASS_7.3-60.0.1                 kSamples_1.2-10                
[217] progress_1.2.3                  uuid_1.2-0                     
[219] BiocParallel_1.36.0             insight_0.19.11                
[221] R6_2.5.1                        rstatix_0.7.2                  
[223] fastmap_1.1.1                   multcomp_1.4-25                
[225] vipor_0.4.7                     rsvd_1.0.5                     
[227] ggstatsplot_0.12.3              nnet_7.3-19                    
[229] gtable_0.3.5                    KernSmooth_2.23-22             
[231] latticeExtra_0.6-30             deldir_2.0-4                   
[233] htmltools_0.5.8.1               RcppParallel_5.1.7             
[235] bit64_4.0.5                     lifecycle_1.0.4                
[237] processx_3.8.4                  spData_2.3.0                   
[239] nloptr_2.0.3                    callr_3.7.6                    
[241] vctrs_0.6.5                     zeallot_0.1.0                  
[243] ggfun_0.1.4                     sp_2.1-4                       
[245] pillar_1.9.0                    gplots_3.1.3.1                 
[247] prismatic_1.1.2                 ShortRead_1.60.0               
[249] locfit_1.5-9.9                  jsonlite_1.8.8                 
[251] expm_0.999-9                    GetoptLong_1.0.5